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Abstract Many reconstruction problems in signal processing require solution of a 
certain kind of nonlinear systems of algebraic equations, which we call Prony systems. 
We study these systems from a general perspective, addressing questions of global 
solvability and stable inversion. Of special interest are the so-called "near-singular" 
situations, such as a collision of two closely spaced nodes. 

We also discuss the problem of reconstructing piecewise-smooth functions from their 
Fourier coefficients, which is easily reduced by a well-known method of K.Eckhoff to 
solving a particular Prony system. As we show in the paper, it turns out that a modi- 
fication of this highly nonlinear method can reconstruct the jump locations and mag- 
nitudes of such functions, as well as the pointwise values between the jumps, with the 
maximal possible accuracy. 



1 Introduction 

In many applications, it is often required to reconstruct an unknown signal from a small 
number of measurements, utilizing some a-priori knowledge about the signal structure. 
Such problems arose (and continue to arise) in recent years under several names in 
different fields, such as Finite Rate of Innovation, super-resolution, sub-Nyquist sam- 
pling and Algebraic Signal Reconstruction [71IHligill3lll4|[]31[Tgi[^ll5T1[M] . One underly- 
ing connection between these problems is that almost all of them require solution of a 
certain kind of nonlinear systems of algebraic equations, which we call Prony systems. 
We therefore consider the study of this system to be an important topic. In particu- 
lar, questions of solvability, uniqueness, including in near-singular situations, as well 
as stability of reconstruction in the presence of noise turn out to be non-trivial and 
requiring a delicate study of some related algebraic-geometric structures. 
This paper consists of two parts. In the first part, we consider the general Prony system. 
First, we present a necessary and sufficient condition for the system to be solvable. Next, 
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we give simple estimate of the stability of inversion in a "regular" setting. Finally, we 
consider inversion in several "near-singular" situations, and in particular the practically 
important situation of colliding nodes. We show that a reparametrization in the basis 
of divided finite differences turns the problem into a well-posed one in this setting. 

In the second part of the paper, we present our recent solution to a conjecture posed by 
K.Eckhoff in 1995 17J, which asks for an algorithm to reconstruct a piecewise-smooth 
function with unknown discontinuity locations from its first Fourier coefficients. While 
the problem of defeating the Gibbs phenomenon received much attention in the last 
decades (see [T1I21I12II17II18II21II22II24II25II33II35| and references therein), the question of 
attaining maximal possible accuracy of reconstruction remained open. We show how 
the Algebraic Reconstruction approach, and in particular an accurate solution of a 
certain Prony system, provides the required approximation rate. 



2 The Prony problem 

Prony system appears as we try to solve a very simple "algebraic signal reconstruction" 
problem of the following form: assume that the signal F(x) is known to be a linear 
combination of shifted <5-functions: 

d 

F (x) = S ajS (x — Xj) . (1) 

3=1 

We shall use as measurements the polynomial moments: 

m k =m k {F) = I x k F(x)dx. (2) 

After substituting F into the integral defining roj. we get 

d d 
171 k {F) = / x k S ' ajS(x — xj) dx = S ' djXj. 

i=i i=i 

Considering aj and unknowns, we obtain equations 

d 

m k (F) = Yj a i xk i> * = 0,1,--.- (3) 

This infinite set of equations (or its part, for k = 0, 1, . . . , 2d — 1), is called Prony 
system. It can be traced at least to R. de Prony (1795, [3D]) and it is used in a wide 
variety of theoretical and applied fields. See [3J for an extensive bibligoraphy on the 
Prony method. 

In writing Prony system @ we have assumed that all the nodes xi, . . . , are pairwise 
different. However, as a right-hand side fj, = (mo, . . . , rn 2 d-i) of is provided by the 
actual measurements of the signal F, we cannot guarantee a priori, that this condition 
is satisfied for the solution. Moreover, we shall see below that multiple nodes may 
naturally appear in the solution process. In order to incorporate possible collisions of 
the nodes, we consider "confluent Prony systems". 
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Assume that the signal F(x) is a linear combination of shifted 5-functions and their 
derivatives: 

s dj-l 

FW^^c/'^-ij). (4) 

3=1 £=0 

Definition 1 For F (x) as above, the vector D (F) = (di, . . . , d s ) is the multiplicity 
vector of F, s = s (F) is its degree and d = Y^,j=i dj is its order. For avoiding ambiguity 
in these definitions, it is always understood that aj ^.^i ^ for all j = 1, . . . , s. 

For the moments = m^(F) = J x k F(x) dx we now get 

3=1 e=o y ' 
Considering and aj % as unknowns, we obtain a system of equations 

S dj ~ 1 k\ 

22 22 -n: ~£)i a j,t x j~ = m fc' = 0, 1, . . . ,2d- 1, (5) 
j=i e=o 

which is called a confluent Prony system of order d with the multiplicity vector D = 
(di, . . . ,d s ). The original Prony system © is a special case of the confluent one, with 
D being the vector (1, . . . , 1) of the length d. 

The system © arises also in the problem of reconstructing a planar polygon P (or 
even an arbitrary semi-analytic quadrature domain) from its moments 

m k{Xp) = II z k xpdxdy, z = x + iy, 
JJr 2 

where \P is the characteristic function of the domain P C M 2 . This problem is im- 
portant in many areas of science and engineering [23] . The above yields the confluent 
Prony system 

s dj — 1 

m k = Cijk(k - 1) • • • (k - i + l)zj~ z , c i . J e C, zj € C \ {0} . 

j=l i=0 

As we shall see below, if we start with the measurements fi(F) = \i = (mo, • • • , wi2d— l)i 
then a natural setting of the problem of solving the Prony system is the following: 

Problem 1 (Prony problem of order d) Given the measurements 

fi = (m , • • ■,rn 2d _ 1 ) e C 2d 

in the right hand side of ©, find the multiplicity vector D = (d%, . . . ,ds) of order 
r = i dj < d, and find the unknowns Xj and aj £, which solve the corresponding 
confluent Prony system ((5]) with the multiplicity vector D. 

It is extremely important in practice to have a stable method of inversion. Many 
research efforts are devoted to this task (see e.g. [4,11, 15,28,29,32"] and references 
therein). A basic question here is the following. 
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Problem 2 (Noisy Prony problem) Given the noisy measurements 

fl = (m , . . .,m 2 d-x) € C 2d 

and an estimate of the error \rhf. — < e^, solve Problem^ so as to minimize the 
reconstruction error. 

3 Solving the Prony problem 

3.1 Prony mapping 

Let us introduce some notations which will be useful in subsequent treatment. 

Definition 2 For each w = (x\, . . . , xj) G C d , let s = s (w) be the number of distinct 
coordinates tj , j = 1, . . . , s, and denote T (w) = (ri, . . . , t s ). The multiplicity vector 
is D = D (w) = (di, ■ ■ ■ , d s ), where dj is the number of times the value Tj appears in 
{x\, . . . , Xd} ■ The order of the values in T (w) is defined by their order of appearance 
in w. 

Example 1 For w = (3,1,2, 1,0,3,2) we have s = 4, T(w) = (3, 1,2,0) and D (w) = 
(2,2,2,1). 

Remark 1 Note the slight abuse of notations between Definition Q] and Definition [5] 
Note also that the order of D (w) equals to d for all w 6 C d . 

Definition 3 For each w 6 C d , let s = s(w), T (w) = (ri,...,r s ) and D (w) = 
(d\, . . . ,d s ) be as in Definition [21 We denote by V w the vector space of dimension d 
containing the linear combinations 

s dj-1 

of 5-functions and their derivatives at the points of T (w). The "standard basis" of V w 
is given by the distributions 

S jti = S W {x-Tj), 3 = 1 «H; I = Q,...,dj-l. (7) 



Definition 4 The Prony space Vd is the vector bundle over C d , consisting of all the 
pairs 

(w,g) : w G C d , g G V w . 

The topology on Vd is induced by the natural embedding Vd C C d x D, where T> is 
the space of distributions on C with its standard topology. 



Finally, we define the Prony mapping VM which encodes the Prony problem. 
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Definition 5 The Prony mapping VM : Vd — ► C 2d for (w, g) £ is defined as 
follows: 

VM((w,g)) = (mo,...,m 2 d-i) eC 2d , m k = m k (g) = J x k g(x)dx. 

Therefore, a formal solution of the Prony problem is given by the inversion of the Prony 
mapping VM. 

Finally, let us recall an important type of matrices which play a central role in what 
follows. 

Definition 6 Let (xi, . . . , x s ) £ C s and D = (di, . . . , d s ) with d = y^ =1 dj be given. 
The d x d confluent Vandermonde matrix is 



VI, V2,0 
VI, 1 V 2 ,l 



V s ,0 
v s,l 



V = V (xi,di, . . . ,x s , d s ) 

- v l,d-l V 2 ,d-1 •■■ V Sid _! 

where the symbol vj ^ denotes the following 1 x dj row vector 

vj,k = f [x$, kx k ~\ . . . ,k (fc - 1) • •• (k - dj) x k r d > +1 ] . 

The matrix V defines the linear part of the confluent Prony system ©, namely, 

mi 



(8) 



V (xi,di, . . . ,x s ,d s ) 



a Mi-i 



a s,d s -l 



(9) 



3.2 Pade problem and the solvability set 

It can be shown that the solution to Problem Q] is equivalent to solving the well- 
known Pade approximation problem. While this connection is extremely important 
and insightful, we do not provide the details here for the sake of brevity. Let us only 
mention the following result. 

Proposition 1 The tuple 



s, D= (di, . . .,d s ), r = ^dj < d, X = {xj}" j=1 , A = {".,.. J _ , _ : 




is a solution to Problem\T\ with right-hand side 

H = (m , . . . ,m 2d _i) £ C 
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if and only if (mo, . . . ,m 2 d-i) are the first 2d Taylor coefficients at z = oa of the 
rational function 

s dj 2d-l 

2d— 1 



*a*.a w = EE « ~ = E5rr+° 

j'=l*=l ^ ^J/ fc=0 



z 



The function Ro.x,A ( z ) * s the Stieltjes transform of the corresponding signal F (x) 

Ej=i Efco 1 ( x - ^j). j - e - 



#d,x,,4 (2) = / 

J — i 



F(x)dt 



Using this correspondence, it is not difficult to prove the following result (see [10|1. 

Theorem 1 Let the right-hand side (mo, . . . , rn-xd— l) °f Problem\T\ be given. Let M d 
denote the dx (d+ 1) Hankel matrix 



M d 



mo mi m 2 
mi m2 ni3 

Wd-i m d m d+1 



^d 

™d+l 
«l2d-l 



For eac/i e ^ d, denote by M e the ex (e + 1) submatrix of M d formed by the first e 
rows and e + 1 columns, and let M e denote the corresponding square matrix. 

Let r ^ d be the rank of M<j. Then Problem\l\ is solvable if and only if the upper left 
minor \M r \ of M d is non-zero. The solution, if it exists, is unique, up to a permutation 
of the nodes {^j}- The multiplicity vector D = (di,...,d s ), y^ =1 dj = r, of the 
resulting confluent Prony system of order r is the multiplicity vector of the poles of the 
rational function Rd,x,A ( z ), solving the Pade problem in Proposition^ 

As a corollary we get a complete description of the right-hand side data n G C 2d for 
which the Prony problem is solvable (unsolvable) . Define for r = 1, . . . , d sets E r C C 2d 
(respectively, E' r C C 2d ) consisting of /i 6 C 2d for which the rank of M d = r and 
\M r \ ^ (respectively, \M r \ = 0). The set E r is a difference X r = X^t \ X 2 of two 
algebraic sets: X* is defined by vanishing of all the s x s minors of M d , r < s < d, 
while X 2 is defined by vanishing of \M r \. In turn, E' r = X r x \X r 2 , with X,. 1 = X* PlX 2 
and X r 2 defined by vanishing of all the r x r minors of M d . The union X r U E' r consists 
of all [i for which the rank of M d = r, which is E\ \ E 2 . 

Corollary 1 The set X (respectively, E r ) of /i G C 2d for which the Prony problem 
is solvable (respectively, unsolvable) is the union X = VJ d =1 E r (respectively, E 1 = 
U d =1 E' r ). In particular, X' C G C 2d ,det M d = 0}. 



So for a generic right hand side fi we have \M d \ ^ 0, and the Prony problem is solvable. 
On the algebraic hypersurface of fi for which \M d \ = 0, the Prony problem is solvable 
if M d _ 1 / 0, etc. 
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3.3 Stable inversion away from singularities 



Consider Problem [5] at some interior point no G S. By definition, no G S ro for some 
ro < d. Let (wo,go) = VM~ X (/xo). Assume for a moment that the multiplicity vector 
Do = D {go) = {di,...d Sa ), X^=l 4j = r o, nas a non-trivial collision pattern, i.e. 
dj > 1 for at least one j = 1, . . . , so. It means, in turn, that the function Rd ,X,A ( z ) 
has a pole of multiplicity dj. Evidently, there exists an arbitrarily small perturbation 
jl of no for which this multiple pole becomes a cluster of single poles, thereby changing 
the multiplicity vector to some D' ^ Do- While we address this problem in Section 
U via the bases of divided differences, in this section we consider a "restricted" Prony 
problem. 

Definition 7 Let PA4 (wo,<7o) = Mo G XV with D {go) = Do and s {go) = so- Let 
Vd denote the following subbundle of P d of dimension so + ro'- 

P Do = {{w,g)eV d : D{g)=D }- 

The restricted Prony mapping VM* D : Pd„ ~ * C So+r ° is the composition 

VM* Dq =koVM \v Dn , 

where ir : C 2d —t C Sa+r ° is the projection map on the first so + ro coordinates. 

Inverting this VM* Dq represents the solution of the confluent Prony system <(5j) with 
fixed structure Do from the first k = 0, 1, . . . ,so + ro — 1 measurements. 

Theorem 2 ([llj) Let Ho = VM* Da {{w ,go)) G C So+r ° with the unperturbed solu- 
tion go = 2?=1 Sfco 1 a j,l^^ > ( x ~ T j)- ^ n a sma M neighborhood of {wo, go) G Pd , 
the map VM*jj is invertible. Consequently, for small enough e, the restricted Prony 
problem with input data fj* 6 C r o+ s o satisfying — /*o|| e ^ as a unique solution. 
The error in this solution satisfies 

i . | 2 (2\ So + ro 1 

where 5 = min,-^- r, — ([for consistency we take aj^i = in the above formula). 

Proof (outline) The Jacobian of PM* D can be easily computed, and it turns out to 
be equal to the product 

JvM* Da = V (n,di + 1, . . . ,T SQ ,d So + l)diag{£,-} 

where V is the confluent Vandermonde matrix (|8j) on the nodes (n , . . . , r So ) and mul- 
tiplicity vector 

D = (di + l,...,d So + 1), 



s 
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while E is the (dj + 1) x (dj + 1) block 

1 ■■■ 
1 ■ ■ ■ aj.o 

000 ... 

Since fj,Q £ S r , the highest order coefficients a^d —1 are nonzero. Furthermore, since all 
the Tj are distinct, the matrix V is nonsingular. Local invertability follows. To estimate 
the norm of the inverse, use bounds from [BJ. □ 

Let us stress that we are not aware of any general method of inverting VM*d q , i.e. 
solving the restricted confluent Prony problem with the smallest possible number of 
measurements. As we shall see below in Section [S] such a method exists for a very 
special case of a single point, i.e. s = 1. 

4 Prony inversion near singularities 

4.1 Collision singularities and finite differences 

Collision singularities occur in Prony systems as some of the nodes in the signal 
F( x ) = y^y— 1 a.iS(x — Xi) approach one another. This happens for /i near the dis- 
criminant stratum A C C 2d consisting of those (mo, • • • , m^d—l) f° r which some of 
the coordinates {zj} in the solution collide, i.e. the function Rd,x,A { z ) nas multi- 
ple poles (or, nontrivial multiplicity vector D). As we shall see below, typically, as /i 
approaches /io € A, i.e. some of the nodes collide, the corresponding coefficients 
a, tend to infinity. Notice, that all the moments = m^(F) remain bounded. This 
behavior creates serious difficulties in solving "near-colliding" Prony systems, both in 
theoretical and practical settings. Especially demanding problems arise in the presence 
of noise. The problem of improvement of resolution in reconstruction of colliding nodes 
from noisy measurements appears in a wide range of applications. It is usually called 
a "super-resolution problem" and a lot of recent publications are devoted to its inves- 
tigation in various mathematical and applied settings. See [13] and references therein 
for a very partial sample. 

Here we continue our study of collision singularities in Prony systems, started in |36| . 
The full details will be published in |10| . Our approach uses bases of finite differences 
in the Prony space Vd in order to "resolve" the linear part of collision singularities. In 
these bases the coefficients do not blow up any more, as some of the nodes collide. 

Let no £ Ed- Consider the noisy Prony problem in a neighborhood of the exact so- 
lution (wo, go) = VM~ X (/^o). As explained in Section 13.31 if D (wo) is non-trivial, 
then there will always be a multiplicity-destroying perturbation, no matter how small 
a neighborhood. Assume that the node vector w = w (/i) is determined, and consider 
the linear system © for recovering the coefficients a = {oj.^} in the standard basis (J7| 
of V w ^y As /x — > fi , the matrix of this linear system will be V (w (/x)) with collision 
pattern D (w (ji)) ^ Do, and therefore its determinant will generically approach zero. 
This will make the determination of {ij^} ill-conditioned, and in fact some of its com- 
ponents will go to infinity. At the limit, however, the confluent problem is completely 
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well-posed since the matrix V (wq, Dq) is non-singular. The challenge is, therefore, to 
make the solution to depend continuously on fi by a suitable change of basis for Vw 
So, instead of the ill-conditioned system 

fi = V(w(jl),D(w (£)))aQ5) 

we would like to have 

£ = V(£)b(£), (10) 

where the matrix V is nonsingular and depends continuously on fj, in the neighborhood 
of fj, . 

First, we extend the well-known definition of divided finite differences to colliding 
configurations. 

Definition 8 Let w = (xi, . . . , xj) be given. For each m = 1,2, ... ,d, denote w m = 
(si, . . . , x m ). According to Definition^ let s m = s (w m ), T (w m ) = (n, m , . . . ,T Sm , m ) 
and D (w m ) = (di,m, • • • , ds m ,m). Consider the decomposition of the rational function 



R w ,m (*) = Yl — — — rf- 

3=1 'j,ml 



into the sum of elementary fractions 

„(m) 



The m-th finite difference A m (w) is the following element of V w : 

s m dj, m (™) 
j=l £=1 

with the coefficients {wff} defined by (fTT]>. 
We prove the following results in |10| . 

Proposition 2 T/ie finite difference Am (w) is a continuous section of the bundle 
Vd- For w G C d with pairwise distinct coordinates, A m (yS) is the usual divided finite 
difference on the elements of w m . 

Theorem 3 For each w £ C d , the collection 



B(w) = {A m {w)} d m=1 C V w 



forms a basis for V w 



Remark 2 Another possible way to construct a good basis B (w) is to build the matrix 
V in (|10p directly by imitating the confluence process of the Vandermonde matrices 
(as done in 20J), multiplying V by an appropriate "divided difference matrix" F (jl) 
(more precisely, a chain of such matrices derived from the confluence pattern). That 
is, V = VF is the new matrix for the recovery of the linear part in (|10|) . while the 
new coefficient vector is b = F _1 a. Also in this case V — > V (wo, Do) as ft — > /^o- The 
matrix F thus defines the corresponding change of basis from {5j^ (w)} as in © to 
B(w). 



10 



D.Batenkov, Y.Yomdin 



Let us now consider the Prony problem in the basis B [w] in some neighborhood of 
HO £ Ed (thus, the order of the exact solution (wo, So) = VM 1 {ho) is d). Writing 
the unknown g € V w in this basis we have 



9 = ^ An An (w) . 



Theorem 4 ([10J) For /i in a sufficiently small neighborhood of no, the solution 
VM- 1 (jl) = (w(jl), {/?m(M)})i 

expressed in the basis B (w) of finite differences, is provided by continuous algebraic 
functions of jl. 



Proof (outline) For each w in a neighborhood of wo, we obtain the system of equations 

d 



^ /3m j x k A m (w) 



m k , k = 0, 1, 



,d- 1. 



(12) 



In the process of solution, the points {asi, . . . , a;^} are found as the roots of the poly- 
nomial Q (z) which appears in the denominator of Rx,D,A ( z ) = q\z] • ^ ne coemc i en t 
vector q of Q (z) is provided by solving a non-degenerate linear system 



M d q: 



m 2 d-i. 



Therefore, w = w (/i) is given by continuous algebraic functions of /i. By Proposition 
[51 the functions 

v k,m ( w ) = / x A m (w)dx 



are continuous in w. At w = wo, the system (|12|l is non-degenerate by assumption, 
therefore it stays non-degenerate in a small neighborhood of mjq. Thus, the coefficients 
fim (w (/t)) are also continous algebraic functions of ft. □ 



4.2 Prony Inversion near £' and Lower Rank Strata 

The behavior of the inversion of the Prony mapping near the unsolvability stratum £' 
and near the strata where the rank of drops, turns out to be pretty complicated. 
In particular, in the first case at least one of the nodes tends to infinity. In the second 
case, depending on the way the right-hand side fi approaches the lower rank strata, the 
nodes may remain bounded, or some of them may tend to infinity. In this section we 
provide one initial result in this direction, as well as some examples. A comprehensive 
description of the inversion of the Prony mapping near £' and near the lower rank 
strata is important both in theoretical study and in applications of Prony-like systems, 
and we plan to provide further results in this direction separately. 



Geometry of Prony systems and reconstruction of piecewise-smooth functions 



11 



Theorem 5 As the right-hand side fj, £ C 2d \ £' approaches a finite point no £ S' , 
at least one of the nodes x%, . . . , x^ in the solution tends to infinity. 

Proof By assumptions, the components mo, . . . , m 2 d-i of the right-hand side /i = 
(mo, • . • ,m 2t ;_i) £ C 2d remain bounded as /i — > no- By Theorem |4j the finite differ- 
ences coordinates of the solution "P.M -1 (/i) remain bounded as well. Now, if all the 
nodes are also bounded, by compactness we conclude that VM (ft) — >" w ^d- By 
continuity in the distribution space (Proposition [2]) we have VM{uj) = no- Hence the 
Prony problem with the right-hand side no has a solution u> £ Vd, in contradiction 
with the assumption that /no € £ . □ 

As it was shown above, for a given n £ E (say, with pairwise different nodes) the rank 
of the matrix is equal to the number of the nodes in the solution for which the 
corresponding 5-function enters with a non-zero coefficients. So n approaches a certain 
fio belonging to a stratum of a lower rank of if and only if some of the coefficients 
aj in the solution tend to zero. We do not analyze all the possible scenarios of such a 
degeneration, noticing just that if 6 ^ , i- e -, the Prony problem is unsolvable for 
^Oi then Theorem [S] remains true, with essentially the same proof. So at least one of 
the nodes, say, ay, escapes to infinity. Moreover, one can show that ajX 2d ~ 1 cannot 
tend to zero - otherwise the remaining linear combination of 5-functions would provide 
a solution for fiQ. 

If no S E, i.e., the Prony problem is solvable for no, all the nodes may remain bounded, 
or some xj may escape to infinity, but in such a way that ajX 2d ~ 1 tends to zero. 

5 Resolution of Eckhoff 's problem 

Consider the problem of reconstructing an integrable function / : [— 7r, 7r] — > R from a 
finite number of its Fourier coefficients 

c fc (/)= f ^- f /(t)e- ,fct di, fc = 0,f,...M. 

^ J — n 

It is well-known that for periodic smooth functions, the truncated Fourier series 

M 

Sm(/)= £<*(/K te 

|fc|=0 

converges to / very fast, subsequently making Fourier analysis very attractive in a vast 
number of applications. By the classical Jackson's and Lebesgue's theorems |27| . if / has 
d continuous derivatives in [— 7T, 7r] (including at the endpoints) and (x) £ Lip (R), 
then 

max \f(x) -5m (/) (*) | < C{R,d) M" d_1 lnM. (13) 

— 7T<X<7T 

Yet many realistic phenomena exhibit discontinuities, in which case the unknown func- 
tion / is only piecewise-smooth. As a result, the trigonometric polynomial $m (/) no 
longer provides a good approximation to / due to the slow convergence of the Fourier 
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series (one of the manifestations of this fact is commonly known as the "Gibbs phe- 
nomenon"). It has very serious implications, for example when using spectral methods 
to calculate solutions of PDEs with shocks. Therefore an important question arises: 
"Can such piecewise-smooth functions be reconstructed from their Fourier measure- 
ments, with accuracy which is comparable to the 'classical' one (1 13 II "? 

It has long been known that the key problem for Fourier series acceleration is the detec- 
tion of the shock locations. Applying elementary considerations we have the following 
fact [5]. 

Proposition 3 Let f be piecewise d-smooth. Then no deterministic algorithm can re- 
store the locations of the discontinuities from the first M Fourier coefficients with ac- 
curacy which is asymptotically higher than M~ d ~ 2 . 

Let us first briefly describe what has become known as the Eckhoff 's method for this 
problem [TollTTirTS] . 

Assume that / has K > jump discontinuities {^/j-.^ (they can be located also at 

±7r, but not necessarily so). Furthermore, we assume that / £ C d in every segment 
l, and we denote the associated jump magnitudes at £j by 

^/=/ ( %+)-/ ( %7)- 

We write the piecewise smooth / as the sum / = F + F, where F(x) is smooth 
and periodic and F(x) is a piecewise polynomial of degree d, uniquely determined by 
{^j } ' { a £J } sucn that it "absorbs" all the discontinuities of / and its first d derivatives. 
This idea is very old and goes back at least to A.N.Krylov (|26|). Eckhoff derives the 
following explicit representation for F(x): 

K d 

j=ie=o (14) 

v * to fc) = - (^njr ( £ ir) ^ ^ + 2 " 

where V n (x; £,•) is understood to be periodically extended to K and B n (x) is the n-th 
Bernoulli polynomial. Elementary integration by parts gives the following formula. 

Proposition 4 Let F(x) be given by (|14[) . For definiteness, let us assume that cq(F) = 
f* n $(x)dx = 0. Then 

K d 

c ^) = ^I>~^I> fc )^~ la ^' ^ = 1,2,.... (15) 
j=i e=o 

Eckhoff observed that if F is sufficiently smooth, then the contribution of Cj : (F) to 
c k{f) is negligible for large k, and therefore one can hope to reconstruct the unknown 
parameters j £j , ag j } from the perturbed equations (| 15p , where the left-hand side reads 
c k (/) ~ c k (^) an d k ^ 1. His proposed method was to construct from the known 
values 

{c fc (/)} k = M - (d + 1) K + 1, M - (d + 1) K + 2, . . . , M 
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an algebraic equation satisfied by the jump points {£1, . . . , £rt}, and solve this equation 
numerically. Based on some explicit computations for d = 1,2; K = 1 and large 
number of numerical experiments, he conjectured that his method would reconstruct 
the jump locations with accuracy M~ . 

We consider the following generalized formulation (without referring to a specific 
method). 

Conjecture 1 (Eckhoff's conjecture) The jump locations of a piecewises-smooth C d 
function can be reconstructed from its first M Fourier coefficients with asymptotic 
accuracy M~ d ~ 2 . 

In [12] we proposed a reconstruction method (see Algorithm [T] below) which is based 
on the original Eckhoff's procedure. 



Algorithm 1 Half-order algorithm, [12) . 



Let / S PC (d, K), and assume that / = + if" where is the piecewise polynomial 
absorbing all discontinuities of /, and 9 6 C d . Assume in addition the following a-priori 
bounds: 



2. Upper bound on jump magnitudes, \ai t j | < A < go. 

3. Lower bound on the value of the lowest-order jump, \ao j | > -B > 0. 

4. Upper bound on the size of the Fourier coefficients of |cj. (<]/)\ < R ■ k~ d ~ 2 . 

Let us be given the first 3M Fourier coefficients of / for M > M (d, K, J, A, B, R) (a quantity 
which is computable). The reconstruction is as follows. 

1. Obtain first-order approximations to the jump locations {£i, . . . , §k} by Prony's method 
(Eckhoff's method of order 0). 

2. Localize each discontinuity £j by calculating the first M Fourier coefficients of the function 
fj=f- hj where hj is a C°° bump function satisfying 

(a) hj = on the complement of [£j — J, £j + J] ; 

(b) hj = lon [fc-f & + £]. 

3. Fix the reconstruction order d± < \^ \ . For each j = 1,2, . . . , K, recover the parameters 
{?j i a 0,j , ■ ■ • , a di .j } from the d\ + 2 equatic 



tions 



i=o v ' 

by Eckhoff's method for one jump (in this case we get a single polynomial equation 
{pm te) = °} of degree di). 

4. From the previous steps we obtained approximate values for the parameters and 
| a; j | . The final approximation is taken to be 

{K di — ^ K di 

3=1 1=0 ) i=l (=0 



We have also shown that this method achieves the following accuracy. 
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Theorem 6 ([12j) Let / 6 PC (d, K) and let f be the approximation of order d\ < 
I 5 J computed by Algorithm\l\ 77ier0 



< d (d, di , K, J, A, B, R) ■ M~ dl 



< C 2 {d,d!,K,J,A,B,R) • M 



l-dx-l 



(16) 



f(x)-f(x) <C 3 (d,di,K,J,A,B,R)-M- 



The non-trivial part of the proof of this result was to analyze in detail the polynomial 
equation p (£j ) = in step 3 of Algorithm [1] It turned out that additional orders of 
smoothness (namely, between d\ and d) produced an error term which, when substi- 
tuted into the polynomial p, resulted in unexpected cancellations due to which the root 
£j was perturbed only by O (K'I~ dl ^. This phenomenon was first noticed by Eckhoff 
himself in 17 for d = 1, but at the time its full significance was not realized. 

An important property of Algorithm [T] is that its final asymptotic convergence order 
essentially depends on the accuracy of step 3. It is sufficient therefore to replace this 



step with another method which achieves full accuracy (i.e. 



M 



-d-2 



) in order to 



obtain the overall reconstruction with full accuracy. It turns out that taking instead of 
consecutive Fourier samples 



k = M — d — 1,M — d. 



,M 



the "decimated" section 



k = N, 2N, . . . , (d + 2) N; 



N 



del' 



M 



(d+2) 



provides this accuracy. 

For the full details, see [5]. Here let us outline our method of proof. 

Denote the single jump point by £ £ [— 7r, it], and let lo = e~^The purpose is to recover 
the jump point lo and the jump magnitudes {ao, . . . , a^} from the noisy measurements 



c fc (/) = ^V— Txr+^fc, k = N,2N,...,(d + 2)N, \e k \ < R ■ k^ 2 . (17) 



Again, we multiply dTTJ by (2n) (ik) d+1 . Denote ctj = i d+1 j a d _j. We get 

d 

m k d =2Tv(ik) d+1 ci : =uj k ^2a j k :i +S k , k = N, 2N, . . . , (d + 2) N, \S k \ < R-k' 1 . 



(18) 



The last (pointwise) bound holds on "jump-free" regions. 
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Definition 9 Let 



„d I. \ def 
PN ~ 



m(j +1)N u 



d+i-j 



Proposition 5 The point u = lu N is a root of p°^ («) 



Proposition 6 The vector of exact magnitudes {ay} satisfies 



m 2 N^ 



m (d+l)A' aJ 



-JV 
-2N 



-(d+l)N 



N 
2N 



N 
(27V) 2 



N 
(2N)° 



1 (d + 1) N {{d + 1) N) . . . {{d + 1) N) c 





~ao~ 








-Ud- 



(19) 



The procedure for recovery of the {cto, ...,a^,u} is presented in Algorithm [2] below, 
while the method for full recovery of the function is outlined in Algorithm [3] below. 



Algorithm 2 Recovery of single jump parameters 



Let there be given the first N 3> 1 Fourier coefficients of the function fj , and assume that the 
jump position £ is already known with accuracy o I N^ 1 ) . 



1. Construct the polynomial 

In («) 



d+l ^ 
3=0 



d+l-j 



from the given noisy measurements mjv, "*2JVi ■ ■ • i m (d+2)N 1181 . 

2. Find the root z which is closest to the unit circle (in fact any root will do). 

3. Take u = \f~z. Note that in general there are N possible values on the unit circle (see 
Remark [3} , but since we already know the approximate location of u> the correct value 
can be chosen consistently. 

4. Recover £ = — argu. 

5. To recover the magnitudes, solve the linear system (19) . 



Remark 3 To see that there are N possible solutions, notice that one recovers e 
which is satisfied by any £ of the form £ = & + ^S"n, n € Z and not just £ = 4r. 



Algorithm 3 Full accuracy Fourier approximation 

Let / E PC (d, K), and assume that / = + where is the piecewise polynomial 
absorbing all discontinuities of /, and £ C . Assume the a-priori bounds as in Algorithm [l] 

1. Using Algorithm ^ obtain approximate values of the jumps up to accuracy 

and the Fourier coefficients of the functions fj. 

2. Use Algorithm P3] to further improve the accuracy of reconstruction. 
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We have shown that indeed full accuracy is acheived. 



Theorem 7 (|5j) Algorithm^ recovers the parameters of a single jump from the given 
noisy measurements ()18[) with the following accuracy: 

< C 4 ^N- d - 2 , 

JD 



<C 5 R(l+^)N- j -\ i = 0,1, 



The main idea of the proof is to analyze the perturbation of the polynomial pfj (u) by 
qff using Rouche's theorem. 

After making the substitution N = [ fg^pgj j i we obtain as an immediate consequence 
of Theorem [7| the resolution of Conjecture [T] 

Theorem 8 Let f € PC (d, K) and let f be the approximation of order d computed 
by Algorithm^ Then 



6-6 



<C 6 (d,K,J,A,B,R)-M 



-d-2 



\ai,j - H,j | < C 7 (d, K, J, A, B, R) ■ M 
J{x)-f(x) <C 8 (d,K,J,A,B,R) ■ M 



l-d-l 



1 = 0,1,. ..,d; 



(20) 



-d-l 



Note that the system (|17[) is a certain variant of the confluent Prony system ([5} for 
just one node. Therefore, Algorithm O can be regarded as a concrete solution method 
for this particular case. 
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